#############################################################################
# Replication code                                                          #
# The Role of Networks in Mobilization for Ethnic Minority Interest Parties #
# By: Simon Otjes                                                           #
# Date: 30/09/2024                                                          #
# Version: 1.0                                                              #
#############################################################################

# Hello! This code comes in five parts
# Part 1: Reading in packages, functions and data
# Part 2: Data preparation
# Part 3: Descriptives
# Part 4: Regression analyses
# Part 5: Visualization

###################################################
# Part 1: Reading in packages, functions and data #
###################################################

remove(list=ls())

###############################
# Part 1.1: overview function #
###############################

overview<-function(vector) {

	return<-c(round(mean(vector,na.rm=TRUE),2),
	round(median(vector,na.rm=TRUE),2),
	round(sd(vector,na.rm=TRUE),2),
	round(min(vector,na.rm=TRUE),2),
	round(max(vector,na.rm=TRUE),2),
	length(vector)-sum(is.na(vector)))
	names(return)<-c("mean","median","sd","min","max","n")
	return(return)
}

##############################
# Part 1.2: read in packages #
##############################

library(foreign)
library(mokken)
library(texreg)
library(Hmisc)
library(margins)

##########################
# Part 1.3: read in data #
##########################

data<-read.spss("~/Dropbox/Data Schaaf Otjes Spierings/replication file/DEMES2021 v1.0.sav")

############################
# Part 2: Data preparation #
############################

################################
# Part 2.1: Dependent Variable #
################################

# Propensity to vote DENK, missing values are greater than 10.

PTV<-as.numeric(data$S101)
PTV[PTV>10]<-NA

DKDENK<-(data$S101=="I don’t know this party")+(data$S101=="DK: Don’t know")

#######################
# Part 2.2: Ethnicity #
#######################

# Four Ethnicity variables: Moroccan-Dutch, Turkish-Dutch, Dutch-Suriname, Dutch-Other Mena and combined Moroccan-Dutch, -Turkish and Other MENA, make sure that true missings are missing

com<-as.numeric(data$N88=="Morocco")
com[data$N88=="NA: Item empty"]<-NA

cot<-as.numeric(data$N88=="Turkey")
cot[data$N88=="NA: Item empty"]<-NA

comena<-as.numeric(data$N88=="Other MENA")
comena[data$N88=="NA: Item empty"]<-NA

comt<-as.logical(com+cot+comena)

#########################
# Part 2.3: Religiosity #
#########################

# religious services, revert to ensure higher values are higher attendance, missings are values above 5

RelSv<-as.numeric(data$V352)
RelSv[RelSv>5]<-NA
RelSv<-6-RelSv

# religious membership: 1 = mentioned; 0 = not mentioned, missings are values above 2

RelMm<-as.numeric(data$V275)-1
RelMm[RelMm>1]<-NA

# religion: islam = 1; other = 0, remove missings

RelDn<-(data$M37=="Islam: Ahmadiyya")+(data$M37=="Islam: Alevism")+(data$M37=="Islam: Salafism")+(data$M37=="Islam: Shi'ism")+(data$M37=="Islam: Sunni")+(data$M37=="Islam: Sufisme")+(data$M37=="Islam: General")+(data$M37=="Islam: Other")
RelDn[data$M37=="DK: Don’t know"]<-NA
RelDn[data$M37=="WS: Won’t say"]<-NA
RelDn[data$M37=="NA: Item empty"]<-NA

#####################
# Part 2.4: Network #
#####################

# network among people without migration background, missings are values above 8, invert so higher values are more contact

netNL<-as.numeric(data$M33)
netNL[netNL>8]<-NA
netNL<-9-netNL

# network among people from country of origin, missings are values above 8, invert so higher values are more contact

netCO<-as.numeric(data$M34)
netCO[netCO>8]<-NA
netCO<-9-netCO

# Network differential

netCN<-netCO-netNL

#partner voted for Ethnic Minority Interest Party, missings are values above 28 

IMEP_ptn<-(data$N130=="NIDA")+(data$N130=="DENK")+(data$N130=="BIJ1")
IMEP_ptn[(data$N130=="INAP: Inapplicable due to routing")|(as.numeric(data$N130)>28)]<-NA
IMEP_ptn[(data$N129=="No")]<-0
IMEP_ptn[(data$N129=="No partner")]<-0

#parents voted for Ethnic Minority Interest Party, missings are values above 3

IMEP_par<-data$N131_8=="Mentioned"
IMEP_par[as.numeric(data$N131_8)>3]<-NA

#children voted for Ethnic Minority Interest Party, missings are values above 3 

IMEP_kid<-data$N132_8=="Mentioned"
IMEP_kid[as.numeric(data$N132_8)>3]<-NA

#family voted for Ethnic Minority Interest Party, missings are values above 3

IMEP_fam<-data$N133_8=="Mentioned"
IMEP_fam[as.numeric(data$N133_8)>3]<-NA

#friends voted for Ethnic Minority Interest Party, missings are values above 3 

IMEP_fri<-data$N134_8=="Mentioned"
IMEP_fri[as.numeric(data$N134_8)>3]<-NA

#colleagues voted for Ethnic Minority Interest Party, missings are values above 3 

IMEP_col<-data$N135_8=="Mentioned"
IMEP_col[as.numeric(data$N135_8)>3]<-NA

#neighbours voted for Ethnic Minority Interest Party, missings are values above 3

IMEP_nei<-data$N136_8=="Mentioned"
IMEP_nei[as.numeric(data$N136_8)>3]<-NA

#scale of the social environment voting for an Ethnic Minority Interest Party 

IMEPNET<-cbind(IMEP_ptn,IMEP_par,IMEP_kid,IMEP_fam,IMEP_fri,IMEP_col,IMEP_nei)

#test scalability (H > 0.3 is acceptable)
#footnote 8 on p.19 and footnote a on p.9 of the appendix
coefH(IMEPNET[rowSums(is.na(IMEPNET))==0,])

# make scale

IMEP_net<-IMEP_ptn+IMEP_par+IMEP_kid+IMEP_fam+IMEP_fri+IMEP_col+IMEP_nei

#################################
# Part 2.7: Use of Social Media #
#################################

# use of social media, missings are values above 5

SocMe<-as.numeric(data$S016)
SocMe[SocMe>5]<-NA

#################################
# Part 2.8: Use of Social Media #
#################################

# left-right distance
# first, position of DENK, missing values are higher than 10.

LiReDENK<-as.numeric(data$V131)
LiReDENK[LiReDENK>10]<-NA

# second, position of self, missing values are higher than 10.

#LiReAll<-mean(LiReDENK,na.rm=TRUE)

LiReSelf<-as.numeric(data$V133)
LiReSelf[LiReSelf>10]<-NA

# absolute distance between self and denk.

LiReDist<-abs(LiReDENK-LiReSelf)

######################
# Part 2.9: Controls #
######################

# Gender: a binary variable 1 = male; 0 = female and non-binary, make sure that true missings are missing

gen<-as.numeric(data$V010=="Male")
gen[data$V010=="NA: Item empty"]<-NA
gen[data$V010=="INAP: Inapplicable due to routing"]<-NA
gen[data$V010=="IC: Did not finish questionnaire"]<-NA

# Age: a binary variable <= 30 = 1; 0 = older, make sure that true missings are missing

age<-(data$V013=="18-20 years")+(data$V013=="21-25 years")+(data$V013=="26-30 years")

age[data$V013=="NA: Item empty"]<-NA

# Education: a binary variable for WO and HBO, make sure that true missings are missing

edu<-(data$V368=="University Master")+(data$V368=="University Bachelor")+(data$V368=="Tertiary higher vocational (i.e., HBO, HTS, HEAO, Kweekschool, Sociale of Pedagogische Academie)")
edu[data$V368=="Not applicable"]<-NA
edu[data$V368=="DK: Don’t know"]<-NA
edu[data$V368=="WS: Won’t say"]<-NA
edu[data$V368=="IC: Did not finish questionnaire"]<-NA
edu[data$V368=="NA: Item empty"]<-NA

############################
# Part 2.10: Make data set #
############################

# make data set

dataset<-as.data.frame(cbind(gen,age,edu,com,cot,comena,comt,netNL,netCO,netCN,IMEP_ptn,IMEP_par,IMEP_kid,IMEP_fam,IMEP_fri,IMEP_col,IMEP_nei,IMEP_net,RelDn,RelSv,RelMm,SocMe,LiReDist,PTV))

########################
# Part 3: Descriptives #
########################

############
# Table A4 #
############

for (i in 1:dim(dataset)[2]){
	print(colnames(dataset)[i])
	print(overview(dataset[,i]))
}

#################################
# Statement about N on pp.17-18 #
#################################

# "Out of the 765 respondents, 361 did not participate in the post-election survey."

length(data$S016)

xtabs(~data$S016=="PA: Did not participate in post-election survey")[2]

##############
# Footnote 5 #
##############

# share of don't know DENK (entire sample)

(xtabs(~DKDENK)/sum(xtabs(~DKDENK)))[2]

# share of don't know DENK (Moroccan-Dutch)

(xtabs(~DKDENK[com==1])/sum(xtabs(~DKDENK[com==1])))[2]

# share of don't know DENK (Turkish-Dutch)

(xtabs(~DKDENK[cot==1])/sum(xtabs(~DKDENK[cot==1])))[2]

# share of don't know DENK (MENA-Dutch)

(xtabs(~DKDENK[comena==1])/sum(xtabs(~DKDENK[comena==1])))[2]

#######################
# Part 4: Regressions #
#######################

#####################################
# Table 1: Regressions in the paper #
#####################################

model01<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+netCN)
model02<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+netCN+SocMe+IMEP_net)
model03<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+netCN+SocMe+SocMe+IMEP_ptn+IMEP_par+IMEP_fri)
model04<-lm(PTV~gen+age+edu+com+comena+RelDn+RelSv*cot+netCN)
model05<-lm(PTV~gen+age+edu+com+comena+RelDn+RelSv*cot+netCN+SocMe+IMEP_net)
model06<-lm(PTV~gen+age+edu+com+comena+RelDn+RelSv*cot+netCN+SocMe+SocMe+IMEP_ptn+IMEP_par+IMEP_fri)

stargazer01<-htmlreg(list(model01,model02,model03,model04,model05,model06), stars = c(0.01, 0.05, 0.1))

######################################
# Table A5: Regressions in the paper #
######################################

# models with one explanatory variables and controls

modelA01<-lm(PTV~com+cot+comena)
modelA02<-lm(PTV~gen+age+edu+com+cot+comena+RelDn)
modelA03<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv)
modelA04<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelMm)
modelA05<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+RelMm)
modelA06<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+SocMe)
modelA07<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+netNL+netCO)
modelA08<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+netCN)
modelA09<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+IMEP_net)
modelA10<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+IMEP_ptn+IMEP_par+IMEP_kid+IMEP_fam+IMEP_fri+IMEP_col+IMEP_nei)
modelA11<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+IMEP_ptn+IMEP_par+IMEP_fri)

stargazerA05<-htmlreg(list(modelA01,modelA02,modelA03,modelA04,modelA05,modelA06,modelA07,modelA08,modelA09,modelA10,modelA11), stars = c(0.01, 0.05, 0.1))

##############################################################
# Table A6: Regression models that split up network variable #
##############################################################

modelA12<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+netNL+netCO)
modelA13<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelSv+netNL+netCO+SocMe+IMEP_net)
modelA14<-lm(PTV~gen+age+edu+com+comena+RelDn+cot*RelSv+netNL+netCO)
modelA15<-lm(PTV~gen+age+edu+com+comena+RelDn+cot*RelSv+netNL+netCO+SocMe+IMEP_net)

stargazerA06<-htmlreg(list(modelA12,modelA13,modelA14,modelA15), stars = c(0.01, 0.05, 0.1))

#################################################################
# Table A7: Regression models with multiple ethnic interactions #
#################################################################

modelA16<-lm(PTV~gen+age+edu+(com+cot+comena)*RelSv+RelDn+netCN)
modelA17<-lm(PTV~gen+age+edu+(com+cot+comena)*RelSv+RelDn+netCN+SocMe+IMEP_net)
modelA18<-lm(PTV~gen+age+edu+(com+cot+comena)*RelSv+RelDn+netCN+SocMe+IMEP_ptn+IMEP_par+IMEP_fri)

stargazerA07<-htmlreg(list(modelA16,modelA17,modelA18), stars = c(0.01, 0.05, 0.1))

########################################################
# Table A8: Regression models with left-right distance #
########################################################

modelA19<-lm(PTV~gen+age+edu+com+cot+comena+RelSv*cot+RelDn+netCN+LiReDist)
modelA20<-lm(PTV~gen+age+edu+com+cot+comena+RelSv*cot+RelDn+netCN+LiReDist+SocMe+IMEP_net)
modelA21<-lm(PTV~gen+age+edu+com+cot+comena+RelSv*cot+RelDn+netCN+LiReDist+SocMe+IMEP_ptn+IMEP_par+IMEP_fri)

stargazerA08<-htmlreg(list(modelA19,modelA20,modelA21), stars = c(0.01, 0.05, 0.1))

####################################################
# Table A9: Regression models for MENA respondents #
####################################################

modelA22<-lm(PTV~gen+age+edu+com+RelSv*cot+RelDn+netCN,data=dataset[comt,])

stargazerA09<-htmlreg(list(modelA22), stars = c(0.01, 0.05, 0.1))

#######################################################################
# Table A10: Regression models with different religions interactions  #
#######################################################################

modelA23<-lm(PTV~gen+age+edu+com+cot+comena+RelDn*(RelSv))
modelA24<-lm(PTV~gen+age+edu+com+cot+comena+RelDn*(RelMm))
modelA25<-lm(PTV~gen+age+edu+com+comena+RelDn+cot*(RelSv))
modelA26<-lm(PTV~gen+age+edu+com+comena+RelDn+cot*(RelMm))

stargazerA10<-htmlreg(list(modelA23,modelA24,modelA25,modelA26), stars = c(0.01, 0.05, 0.1))

####################################################################
# Table A11: Regression models with interactions for social media  #
####################################################################

modelA27<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+SocMe*(netNL+netCO+IMEP_net))
modelA28<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+SocMe*(netCN+IMEP_net))
modelA29<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+SocMe*(netCN+IMEP_ptn+IMEP_par+IMEP_fri))

stargazerA11<-htmlreg(list(modelA27,modelA28,modelA29), stars = c(0.01, 0.05, 0.1))

#####################
# Save regressions  #
#####################

write(stargazer01,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table 01.html")
write(stargazerA05,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A05.html")
write(stargazerA06,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A06.html")
write(stargazerA07,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A07.html")
write(stargazerA08,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A08.html")
write(stargazerA09,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A09.html")
write(stargazerA10,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A10.html")
write(stargazerA11,"~/Dropbox/Data Schaaf Otjes Spierings/replication file/Tables/Table A11.html")

###################
# Part 5: Figures #
###################

########################
# Effect sizes on p.23 #
########################

# "among those who attend religious services once a week or more, most have a PTV nearly one point higher than those who never do"

pred_dat21 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelSv = c(1,5))

pred21 <- predict(model01,newdata=pred_dat21, interval = "confidence", level = 0.95)

round(abs(pred21[1,1]-pred21[2,1]),1)

# "Model A4 shows that respondents who are members of a religious organization have a one-and-a-half points higher PTV for DENK than those who are not."

pred_dat22 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelMm = c(0,1))

modelA04<-lm(PTV~gen+age+edu+com+cot+comena+RelDn+RelMm)

pred22 <- predict(modelA04,newdata=pred_dat22, interval = "confidence", level = 0.95)

round(abs(pred22[1,1]-pred22[2,1]),1)

#footnote 14: "the explanatory power of the model with the interaction effect of attendance for being Turkish is slightly higher than for being Muslim (Table A10)"

summary(modelA25)$r.squared-summary(modelA23)$r.squared
summary(modelA26)$r.squared-summary(modelA24)$r.squared

# "For them, the increase is three points on the PTV for DENK, compared to a non-significant one-point effect found for other respondents."

pred_dat23 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 1,
		RelSv = c(1,5))

pred23 <- predict(model04,newdata=pred_dat23, interval = "confidence", level = 0.95)

round(abs(pred23[1,1]-pred23[2,1]),1)

pred_dat24 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelSv = c(1,5))

pred24 <- predict(model04,newdata=pred_dat24, interval = "confidence", level = 0.95)

round(abs(pred24[1,1]-pred24[2,1]),1)

#significant?
round(abs(pred24[1,1]-pred24[2,1]),1)>((pred24[1,2]-pred24[1,3])^2+(pred24[2,2]-pred24[2,3])^2)^(1/2)

#"The effect is clearly positive, as expected, with a maximum impact of more than three points on the PTV."

pred_dat25 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = c(-7,7),
		cot = 0,
		RelSv = 1.857)

pred25 <- predict(model04,newdata=pred_dat25, interval = "confidence", level = 0.95)

round(abs(pred25[1,1]-pred25[2,1]),1)

#########################################
# Figure 1: Four ethnicity interactions #
#########################################

#########################
# Figure 1a: set values #
#########################

pred_dat01 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 1,
		RelSv = c(1:5))

####################################
# Figure 1a: calculate predictions #
####################################

pred01 <- predict(model04,newdata=pred_dat01, interval = "confidence", level = 0.95)

#########################
# Figure 1b: set values #
#########################

pred_dat02 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 1,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelSv = c(1:5))

####################################
# Figure 1b: calculate predictions #
####################################

pred02 <- predict(model04,newdata=pred_dat02, interval = "confidence", level = 0.95)

#########################
# Figure 1c: set values #
#########################

pred_dat03 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 1,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelSv = c(1:5))
		
####################################
# Figure 1c: calculate predictions #
####################################		

pred03 <- predict(model04,newdata=pred_dat03, interval = "confidence", level = 0.95)

#########################
# Figure 1d: set values #
#########################

pred_dat04 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = -0.4745,
		cot = 0,
		RelSv = c(1:5))

####################################
# Figure 1d: calculate predictions #
####################################	

pred04 <- predict(model04,newdata=pred_dat04, interval = "confidence", level = 0.95)

# four figures

par( mfrow = c( 2, 2 ) )

###################
# Figure 1a: plot #
###################

bp<-barplot(xtabs(~RelSv),col="white",border="grey",xlab="Religious Attendance",ylab="PTV DENK",main="Turkish-Dutch voters",yaxt="n",ylim=c(0,600))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*600,labels=c(0,2,4,6,8,10))

lines(x=bp,y=pred01[,1]*60)	
lines(x=bp,y=pred01[,2]*60,lty=2)	
lines(x=bp,y=pred01[,3]*60,lty=2)	

###################
# Figure 1b: plot #
###################

bp<-barplot(xtabs(~RelSv),col="white",border="grey",xlab="Religious Attendance",ylab="PTV DENK",main="Moroccan-Dutch voters",yaxt="n",ylim=c(0,600))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*600,labels=c(0,2,4,6,8,10))

lines(x=bp,y=pred02[,1]*60)	
lines(x=bp,y=pred02[,2]*60,lty=2)	
lines(x=bp,y=pred02[,3]*60,lty=2)	

###################
# Figure 1c: plot #
###################

bp<-barplot(xtabs(~RelSv),col="white",border="grey",xlab="Religious Attendance",ylab="PTV DENK",main="MENA-Dutch voters",yaxt="n",ylim=c(0,600))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*600,labels=c(0,2,4,6,8,10))

lines(x=bp,y=pred03[,1]*60)	
lines(x=bp,y=pred03[,2]*60,lty=2)	
lines(x=bp,y=pred03[,3]*60,lty=2)	

###################
# Figure 1d: plot #
###################

bp<-barplot(xtabs(~RelSv),col="white",border="grey",xlab="Religious Attendance",ylab="PTV DENK",main="Other voters",yaxt="n",ylim=c(0,600))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*600,labels=c(0,2,4,6,8,10))

lines(x=bp,y=pred04[,1]*60)	
lines(x=bp,y=pred04[,2]*60,lty=2)	
lines(x=bp,y=pred04[,3]*60,lty=2)	

#######################################
# Figure 2: Embeddedness differential #
#######################################

########################
# Figure 2: set values #
########################

pred_dat05 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN = c(-7:7),
		cot = 0,
		RelSv = 1.857)

###################################
# Figure 2: calculate predictions #
###################################

pred05 <- predict(model04,newdata=pred_dat05, interval = "confidence", level = 0.95)

##################
# Figure 2: plot #
##################

bp<-barplot(xtabs(~netCN),col="white",border="grey",ylab="PTV DENK",xlab="Embeddedness Differential",main="",yaxt="n",ylim=c(0,300))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*300,labels=c(0,2,4,6,8,10))

lines(x=bp,y=pred05[,1]*60)	
lines(x=bp,y=pred05[,2]*60,lty=2)	
lines(x=bp,y=pred05[,3]*60,lty=2)	

##############################
# Figure 3: Social media use #
##############################

########################
# Figure 3: set values #
########################

pred_dat06 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=-0.4745,
		cot = 0,
		RelSv = 1.857,
		SocMe=c(1:5),
		IMEP_net=0.6038)
		
###################################
# Figure 3: calculate predictions #
###################################		

pred06 <- predict(model05,newdata=pred_dat06, interval = "confidence", level = 0.95)

##################
# Figure 3: plot #
##################

bp<-barplot(xtabs(~SocMe),col="white",border="grey",xlab="Politics news on social media",ylab="PTV DENK",main="",yaxt="n",ylim=c(0,250))
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*250,labels=c(0,2,4,6,8,10))
lines(x=bp,y=pred06[,1]*25)	
lines(x=bp,y=pred06[,2]*25,lty=2)	
lines(x=bp,y=pred06[,3]*25,lty=2)	

#################################
# Figure 4: IMEP network voting #
#################################

########################
# Figure 4: set values #
########################

pred_dat07 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=-0.4745,
		cot = 0,
		RelSv = 1.857,
		SocMe=3.602,
		IMEP_net=0:6)

###################################
# Figure 4: calculate predictions #
###################################	

pred07 <- predict(model05,newdata=pred_dat07, interval = "confidence", level = 0.95)

##################
# Figure 4: plot #
##################

bp<-barplot(xtabs(~IMEP_net),col="white",border="grey",xlab="Network EMIP voting",ylab="PTV DENK",main="",yaxt="n",ylim=c(0,250)*1.2)
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*250,labels=c(0,2,4,6,8,10))
lines(x=bp,y=pred07[,1]*25)	
lines(x=bp,y=pred07[,2]*25,lty=2)	
lines(x=bp,y=pred07[,3]*25,lty=2)

#######################################################
# Figure A1: Social Media * IMEP Network interactions #
#######################################################

##############################
# Figure A1: set both values #
##############################

pred_datA1 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=-0.4745,
		cot = 0,
		RelSv = 1.857,
		SocMe=c(1:5),
		IMEP_net=0)

pred_datA2 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=-0.4745,
		cot = 0,
		RelSv = 1.857,
		SocMe=c(1:5),
		IMEP_net=2)

####################################
# Figure A1: calculate predictions #
####################################

predA1 <- predict(modelA28,newdata=pred_datA1, interval = "confidence", level = 0.95)
predA2 <- predict(modelA28,newdata=pred_datA2, interval = "confidence", level = 0.95)

###################
# Figure A1: plot #
###################

bp<-barplot(xtabs(~SocMe),col="white",border="grey",xlab="Politics news on social media",ylab="PTV DENK",main="",yaxt="n",ylim=c(0,250)*1.2)
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*250,labels=c(0,2,4,6,8,10))
lines(x=bp,y=predA1[,1]*25,col="grey")	
lines(x=bp,y=predA1[,2]*25,col="grey",lty=2)	
lines(x=bp,y=predA1[,3]*25,col="grey",lty=2)

lines(x=bp,y=predA2[,1]*25)	
lines(x=bp,y=predA2[,2]*25,lty=2)	
lines(x=bp,y=predA2[,3]*25,lty=2)

##############################
# Figure A2: set both values #
##############################

pred_datA3 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=-4,
		cot = 0,
		RelSv = 1.857,
		SocMe=c(1:5),
		IMEP_net=0.6038)

pred_datA4 <- data.frame(
		gen = 0,
		age = 0,
		edu = 0,
		com = 0,
		comena = 0,
		RelDn = 0,
		netCN=1,
		cot = 0,
		RelSv = 1.857,
		SocMe=c(1:5),
		IMEP_net=0.6038)

####################################
# Figure A2: calculate predictions #
####################################

predA3 <- predict(modelA28,newdata=pred_datA3, interval = "confidence", level = 0.95)
predA4 <- predict(modelA28,newdata=pred_datA4, interval = "confidence", level = 0.95)

###################
# Figure A2: plot #
###################

bp<-barplot(xtabs(~SocMe),col="white",border="grey",xlab="Politics news on social media",ylab="PTV DENK",main="",yaxt="n",ylim=c(0,250)*1.2)
axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1.0)*250,labels=c(0,2,4,6,8,10))
lines(x=bp,y=predA3[,1]*25,col="grey")	
lines(x=bp,y=predA3[,2]*25,col="grey",lty=2)	
lines(x=bp,y=predA3[,3]*25,col="grey",lty=2)

lines(x=bp,y=predA4[,1]*25)	
lines(x=bp,y=predA4[,2]*25,lty=2)	
lines(x=bp,y=predA4[,3]*25,lty=2)

###########
# The End #
###########